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"Method for determining the three-dimensional surface of an object". 

* * * * 

DESCRIPTION 

The present invention refers to a method for determining the three- 
5 dimensional surface of an object and relevant computer program. 

For describing and analysing the geometry of objects that have to be 
produced, modelling methods that use computers are known. In particular, 
when a three-dimensional surface enclosed by a set of points acquired from 
a real object has to be obtained, a technique called level-set that describes an 

10 enclosed surface as a level surface of a volumetric function (or implicit 

function) is normally used. 

The level-set technique provides for the time evolution of the 
volumetric function, according to a suitable equation at the partial 
derivatives, typically belonging to the Hamilton - Jacobi family of equations. 

1 5 The volumetric function describes the distance of each point from the 

enclosed surface and this distance will have positive or negative sign 
according to the point considered whether it is inside or outside the surface. 
Thus the set of points that possess a null distance from the surface S 
represent the surface itself. The technique provides for the evolution of the 

20 enclosed surface, defined by flie Hamilton - Jacobi equation so that it 

surrounds and adheres to the cloud of 3D points that describes the object. 

This technique is used as it is capable of adapting to the various types of 
objects and to the most diverse point acquisition techniques. This technique 
in fact, does not require any additional information about the topology of the 

25 points acquired and is insensitive to the technique used: the density of the 

points can vary in the different zones of the object considered; no 
information is required about the point acquisition sequence or about their 
relative spatial positions; and the number of objects enclosed separated from 
each other is not relevant for the algorithm. 

30 As this technique requires the updating of a volumetric function at each 
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step of evolution of the front, the calculation cost is usually very high, so 
much so that it makes the method of very little practical utility. Methods are 
known that permit the reduction of the calculation cost of this technique, for 
example limiting the updating of the volumetric function in an enclosed 
5 zone surrounding the front in evolution. Despite this the calculation cost 

remains prohibitive. 

Techniques in alternative to the level-set technique are also known. For 
example, methods such as those of the Delaunay tri angulation are known, as 
well as the NURBS (Non Uniform Rational B -Splines) and the HRBF 
10 (Hierarchical Radial Basis Function). Each of these methods has limitations 

and impose a series of constrictions on the set of the points acquired, such as 
for example the calculation complexity depends heavily on the number of 
points, or they need to have a spatial distribution of the points acquired as 
uniform as possible. 

15 In view of the state of the technique described, object of the present 

invention is to provide for a rapid method for determining the three- 
dimensional surface of an object and that does not have the inconveniences 
of the known art. 

In a first aspect the present invention relates to a method for 

20 determining the three-dimensional surface of an object comprising the 

phases of: defining the coordinates of a plurality of points of said object; 
defining a three-dimensional matrix of cells that contains said object to 
which a value can be associated; determining the distance between each 
centre of said cells of said three-dimensional matrix of cells and the closest 

25 point of said plurality of points of said object; setting the value of several 

cells of said three-dimensional matrix of cells at a first preset value; 
determining the value that each cell of said three-dimensional matrix of cells 
assumes, with the exception of said several cells, by means of the following 
formula 



WO 2005/027053 



PCT/EP2004/052154 



3 



F(x i9 t) • v, + w • Y,F(x J9 f) • v, 

where 

3c. represents the coordinates of the centre of the i_th cell, 
F(x.,t) represents the value of the i_th cell at time t, 
v, represents said distance, 
5 w represents a second preset value, and 

j indicates a neighbourhood of cells of the i_th cell; 
determining the module sum of the variations of the value of each cell 
between the time / and the time t+1; repeating said phase of determining the 
value that each cell of said three-dimensional matrix of cells assumes if said 
1 0 sum is higher than a third preset value. 

In its second aspect the present invention relates to a computer program 
comprising a program code that carries out all the phases of the method for 
determining the three-dimensional surface of an object when said program is 
carried out on said computer. 
15 In its third aspect the present invention relates to a computer program 

memorised on a support that can be used by said computer to control 1he 
execution of all the phases of the method for detennining the three- 
dimensional surface of an object 

Thanks to the present invention, by substituting the differential equation 
20 at the partial derivatives normally used in literature (Hamilton- Jacobi) with 

one that is very well suited to describing discontinuous functions such as 
those that describe the interface between two immiscible fluids, that is one 
of the Navier-Stokes equations, derived from the fluid dynamics, the 
calculation cost is considerably reduced, and a much quicker convergence 
25 can be obtained. The surface can be interpreted as the front between two 

viscous fluids that moves until it completely encloses the object. In 
particular, acting on the viscosity value of the fluids used it is also possible 
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to control the level of roughness permitting smooth surfaces to be obtained 
even if the points presented deviations from their correct spatial position. In 
addition, the method proposed permits the introduction of elements and 
strategies that enable an end control of the evolution of the front, based on 
5 principles of simple physics interpretability. For example, an artificial 

turbulence can be introduced in the fluids that facilitates penetration inside 
particularly narrow cavities. Thus it is possible to provide a correct 
determination even of configurations of points that traditionally had been 
negated. Thanks to its flexibility, the Navier-Stokes equation can be used for 

10 determining more volumetric functions that describe different parts of the 

object considered. Objects that are particularly complex or that have high 
resolutions can thus be represented by calculating separately the evolution of 
the volumetric function on the different parts that make them up. The single 
results obtained are then recomposed to obtain the overall volumetric 

1 5 function using a simple mathematic operator. 

The characteristics and advantages of the present invention will appear 
evident from the following detailed description of an embodiment thereof, 
illustrated as non- limiting example in the enclosed drawing, in which Figure 
1 shows a flow diagram of the method for determining the three-dimensional 

20 surface of an obj ect. 

The present invention is based on the Navier- Stokes differential 
equation at the partial derivatives for the conservation of mass in which a 
redefinition of the velocity vector is used according to a criterion described 
hereunder. It describes the motion of two immiscible fluids that converge, 

25 one from inside and the other from outside, towards tihe cloud of points that 

describes the object considered. At the end of the evolution, when both the 
fluids adhere to the points of the object, the surface of separation between 
the two coincides with the surface of the object itself. 

The volumetric function describes the contents of the two fluids in each 

30 elementary cell or voxel, where voxel means the elementary element of the 
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volumetric grid, such as how the pixel is used for a two-dimensional grid. 
To distinguish the two fluids we have assigned opposite sign mass to them. 
Therefore a cell in which the volumetric function is worth -1 contains 
exclusively the internal fluid, while a value of 4-1 indicates that only external 
5 fluid is present 

If between two adjacent cells the volumetric fiinction changes sign it 
implies that the surface of separation passes between them and thus also the 
surface of the object. The exact point of passage is determined by means of 
interpolation operations. 
10 The exact position of the contact surface between the two fluids is thus 

identified by the position of the level-set zero (that is the set of points with 
null distance from the surface). 

The law for the conservation of mass is one of the fundamental 
principles of classic physics and it is independent from the nature of the 
15 fluids and from the forces that act on it; it expresses the principle according 

to which in a fluid dynamic system mass cannot disappear or be created 
unless sinks or sources are present. 

The value F of the volumetric function in each cell indicates the 
quantity of fluid inside it and this value is modified by the forces acting on 
20 the fluid, by its velocity, by the diffusion coming from the adjacent cells and 

from the sources or sinks. 

The evolution of a fluid inside a system can be described by means the 
flow vector G that expresses the quantity of fluid that crosses a surface in a 
given interval of time. 
25 In its classic formulation the law of conservation for mass affirms that 

the variation by time unit of the quantity of matter F inside a volume £2, 

d 



UFda 



a 

is equal to the net contribution of the flow that crosses the external 
surface S of Q, to the contribution Qv of sources (or sinks) of volume inside 
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it, and to the contribution of sources Q s orientated, along its surface S. The 
contributions due to the different sources can thus be written as: 

a s 

The general formulation for the conservation of mass equation can thus 
~\ Fd£l = -|G dS+\Q v da+\Q 5 • dS 

dt a S ft s 

be written as: 

The negative sign in front of the flow integral is due to the feet that the 
area dS is orientated outwards from £2. 

Using the Gauss theorem the previous equation can be rewritten using 
exclusively volume terms: 



dF 
dt 

The flow density vector G can be decomposed in a diffusive 
10 component G D and in a convective component G c , the convective part of 

the flow vector describes the transport phenomenon due to the external 
forces and is defined as product between the velocity v of the flow and the 
quantity of matter transported F: therefore it results G c =vF . 

The diffusive component instead derives from the remixing, present also 
15 in resting fluids, due to the thermal agitation, and is usually proportional, 

according to the law of Fick, to the gradient of F 

G D =yVF where y is the diffusivity constant The time evolution of the 
volumetric function F due to the diffusivity is regulated (as can be verified 
forcing the conservation of mass) by the differential equation at the partial 
20 derivatives of the heat and thus its contribution results mainly connected to 

regularizing and a flattening of the function F. The effect of the diffusive 
contribution is thus taken into consideration applying a Gaussian filtering to 
the entire volumetric function. 

The intensity and the orientation of the velocity v of the fluid is 
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determined in each cell of the matrix. This velocity has been defined 
orientating the flow towards the point of the 3D cloud (that describes the 
object) closest to the cell, and the module of this vector has been selected 
directly proportional to the distance: in this manner both the fluids 
progressively slow down coming closer to the cloud of points and converge 
at the final surface of contact The equation that has thus been used to 
describe the evolution of the two fluids is the following: 



fF(x,0lv(x)|-e 2 ° 2 dx 

dt r, , ^¥ 

j|v(jc)| e 2a dx 

a 

The local variation of the quantity of fluid contained in the i_th cell 
F(x i9 t) between two successive instants therefore results determined by the 

10 convective contribution deriving from each cell F(x 9 t)\v(x)\ weighed with a 

Gaussian function whose standard deviation o describes the degree of 
diffusion of the fluid between the different cells. 

The velocity v (3c) has been defined for each cell proportional to the 
distance between the centre of the cell and the closest point belonging to the 

15 3D cloud. In particular the module of this velocity equals I*,. - /?| a where x { 

is the centre of the ijh cell, while p is the closest point of the cloud of 
points (the point at least distance from the centre of the i_th cell), a is 
instead a parameter that regulates the course of the velocity function. 
Experimentally it has been determined that acceptable values for a are 

20 between 1,5 and 2,5. The greater this value is, the slower the velocity will be 

in proximity of the points while it will be rapid at a great distance. Even 
though this characteristic is very important for a rapid evolution of the 
values of the matrix it can however lead to a slow final convergence. A 
preferred value of the value of a equals 2.0. 

25 Another property of the fluids used in the model is the viscosity, whose 
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value determines the standard deviation o inside the equation. This property 
describes the viscose friction that acts inside the fluid limiting its mobility. 
Considering this property we prevent the fluids that cross the reconstructed 
surface from flowing between the interstices of the various points. This 
5 property is particularly important as through it special arrangement of points 

can be correctly met such as fine blades or particularly pointy objects. 
Thanks to the viscosity, the fluids do not cross the surface which remains 
adherent to the blade or to the point. It is also possible to reduce the curving 
of the entire surface obtaining particularly smooth surfaces and reducing the 

1 0 surface roughness that can also have originated from the error in 

measurement. It must be remembered however that high values can however 
limit the capacity of the algorithm to represent objects with particularly 
jagged surfaces and to enable the flow of the liquid inside concave zones. 
Because of the above-mentioned the viscosity is a parameter whose choice is 

1 5 determined according to the typology of the object considered, by the 

imprecision of the acquisition system used and by the degree of definition 
required for the adhesion of the final surface to the cloud of points. 

Another aspect that has been used by the physical modelling of the 
fluids is the turbulence, which treats the presence of whirlpools that describe 

20 local variations of the velocity function. 

The presence of the turbulence can locally develop high pressures and 
can generate instability in the flow but at the same time can favour the 
penetration of the fluid inside very small openings and prevent stagnation in 
zones with modest velocities. 

25 The presence of vortexes was introduced modulating the field of 

velocity with small oscillations in the different directions around the initial 
configuration. Because of the close correlation between the velocity vector 
in each cell and the vector that points towards the closest point of the 3D 
cloud, the local oscillations of the velocity field can also be interpreted as an 

30 elastic deformation of the cloud of points in the various directions which can 
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facilitate the entrance of the fluid inside the restricted zones. 

The modulation of the velocity function is applied cyclically along the 

three spatial axes multiplying by a factor the terms v of the velocities of the 

cells that lie orientated in parallel to these axes in relation to the i_th cell 
5 under examination and its intensity can be controlled externally, in particular 

at each iteration the velocity is amplified in a given direction and/or the 

velocity in the other two directions is reduced by an equal quantity. 

Obviously greater the expansion coefficients used greater will be the 

turbulence. Values that are too high can however determine an excessive 
10 remixing also between the two fluids on the contact surface causing the 

leakage of the internal fluid or an excessive penetration of the external fluid 

in the object considered. 

The previous equation that was used to describe the evolution of the two 

fluids has been discretized to be able to carry out the operations more easily 
15 by means of a computer and is described by the following equation. 

F(x i9 1 + 1) - F(x i9 1) = £ F(x i9 t) 

J 

where 

J. represents the spatial coordinates (three-dimensional array) of the 
ijhcell, 

F(x i9 t) represents the value of the i_th cell at time t, 
20 v, represents the velocity (or the distance as previously defined), 

w represents a preset value of viscosity, and 

j is a preset value that indicates a neighbourhood of cells of the i_th 

cell. 

In particular, the index summation j considers the contribution of the 
25 closest cells to the point under examination. Preferably, this interaction was 

chosen to be limited exclusively to the first 6 cells near each cell under 



WO 2005/027053 



PCT/EP2004/052154 



10 



examination, that is that share with it an entire face of the cubical cell (the 
cells at distance 1 from the cell under examination). 

We now refer to Figure 1 that shows the flow diagram of the method for 
detennining the three-dimensional surface of an object This flow diagram 
5 describes how a computer program carries out the calculation for 

determining the three-dimensional surface of an object. 

Phase 1 - Data definition. 

In this phase the coordinates of a plurality of points of the object that is 
needed to be represented are defined. The set of the points is described as a 
10 set of vectors with 3 components, containing the spatial coordinates of each 

point No additional information is requested in relation to the points. 

Phase 2 - Definition of cell matrix. 

In this phase the cell matrix (or three-dimensional grid) that defines the 
evolution dominium of the value of F is described. The number of cells is 

15 chosen on the basis of the dimensions of the object that is needed to be 

represented. The number of cells that will be used is determined on the basis 
of the extensions along the three axes of the object Thus keeping a cubical 
form for each elementary cell of the matrix, a grid that contains exactly out 
object is defined. It is also possible to define a form of the cells that is not 

20 cubical, that is a different resolution of the matrix along the different axes 

can be set. This can be very important for objects in which a high resolution 
is necessary only in a dimension such as low relief, building feces, etc.. 

It is also possible to define a rotated matrix in relation to the axes 
according to which the object has been acquired and, in particular, align it 

25 with the main directions of the object itself. Thus it is possible to minimize 

the number of cells used at the same time keeping the same resolution. A 
centre with the appropriate coordinates belongs to each cell, and a value can 
be associated to each cell. 

Phase 3 - Detennining the velocity. 

30 In this phase the field of the velocities that belong to each cell of the 



WO 2005/027053 



PCT/EP2004/052154 



11 



grid is defined. As we saw previously the velocity field depends closely on 
the distance of each centre of the cells from the closest point of the cloud of 
points. Thus, the field of the distances is a three-dimensional matrix with the 
same dimensions as the grid. For each cell the point of the 3D cloud that is 
5 the closest (at least distance) is determined and the distance is calculated as 

|jc, -pf between this point and the centre of the cell. This value is 
memorized in equivalent cell of the matrix of the distances. 

Take note that it is not necessary to know the direction of the velocities 
but only the module. 

10 It is important to note how, the dependence of the calculation 

complexity of this method on the number of points of the 3D cloud is 
limited exclusively to the determination of the velocity field, only once at 
the beginning of the procedure and all the evolutive phase it is completely 
independent from the number of sample points used. 

1 5 Phase 4 — Arrangement of the sources. 

Different arrangements of the sources are possible. A first arrangement 
provides for the external sources to be placed exclusively on the edges of the 
matrix (which will thus always be kept at the value +1 during the entire 
evolution) and all the remaining cells instead to be filled with external fluid 

20 (value equal to - 1). A second arrangement instead provides also for one or 

more internal sources to be provided, in this case the cells of the edge will 
always be kept at +1, the value of the internal cells will be equal to 0, while 
the internal sources will always be kept at -1 . This method permits an even 
more rapid evolution but requires initial information relative to the 

25 arrangement of the additional internal sources that are not easily available. 

Phase 5 — Fixing viscosity and turbulence parameters. 
The viscosity (in the discretized formula represented by w) regulates the 
diffusion of the fluid and thus the weight of the interaction between each 
single cell and its neighbours. The preferred value of w results between 0,1 

30 and 0,9. 
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The parameters relating to the turbulence indicate how the velocities of 
the fluid, in each cell, are altered at each iteration with the purpose of 
simulating the turbulent effect In particular it is possible to set for each of 
the 3 axes the coefficient by which the velocity in this direction will be 
5 multiplied cyclically. The multiplicative factors for the velocity vary 

preferably by a minimum of 0.5 up to a maximum of 1.5. The use of the 
turbulent regime is optional. 

Phase 6 - Starting up the procedure. 

To accelerate the operations, and optionally, provision is made in this 
10 field for another three-dimensional matrix with the dimensions of the matrix 

of the cells in which the information regarding the updating or not of the 
different cells is contained. If a cell can be updated its contents will be 
processed by the evolutive algorithm and thus the contents of liquid inside it 
can vary between two iterations. Cells that cannot be updated are the sources 
15 (that are always kept at the same value, +1 or -1). During the evolution also 

those cells for which the variation of liquid between two iterations has not 
been relevant (below a certain set threshold) can become non-updateable and 
thus these cells are frozen until the evolutive front approaches them. 
Phase 7 - Updating cells. 
20 It is the main phase of the method in which the value of the quantity of 

fluid present in each cell of the three-dimensional matrix of cells is 
determined. 

The use of the Gaussian function inside the integral would provide for 
the analysis at each iteration of the interaction of each cell with all the 
25 matrix of the cells, which would result in being particularly cumbersome in 

computation terms. 

In the case with the discretized formula, it is possible and preferable to 
limit the extension of this interaction exclusively to a neighbourhood of cells 
closer to the i_th cell, that can be for example 6 or 18 or 26 (the cell 
30 considered is at the centre of a matrix 3x3x3), that is those that are 
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respectively in face, edge or vertex contact (that is at distance 1, V2 and V3 
from the cell under examination, in the case of cubic cells) from the i_th 
cell. Preferably the 6 closest cells have been chosen, that is the 6 cells that 
contact the 6 faces of the i_th cell. 
5 Once the cell has been updated, the variation in relation to the previous 

value is assessed and if this value is lower than a threshold set by the user 
the contents of the cell are frozen and the cell is declared not updateable. 
Phase 8 - Visualization 

It is possible, if it interests, to visualize the level-set zero, that is the 
10 surface of separation between the two fluids. The algorithm used is that 

commonly known with the name of Marching Cubes (see for example the 
article "Marching cubes: a high resolution 3D surface construction 
algorithm", by William E. Lorensen and Harvey E. Cline, from Computer 
Graphics, Volume 21, Number 4, July 1987) by means of which it is 
15 possible to obtain an effective triangulation of the level-set zero and a sub- 

pixel resolution thanks to the linear interpolation carried out on adjacent 
cells. The triangulation thus obtained is represented on video, for example, 
by means of an interface program called OpenGL that permits visualization 
in real time and a direct interaction with the user who can interactively 
20 change at will the view point analysing the evolution of the level-set in real 

time. 

Phase 9 - Determining the flow variations. 

At the end of each evolution the overall quantity of fluid that has been 
moved is assessed. 
25 Phase 10 - Evolution cycle. 

If the overall quantity of fluid that has been moved is slower than a 
preset threshold the evolution is interrupted at phase 1 1, otherwise it carries 
out further iteration. 

At phase 1 1, when the calculation is completed, a three-dimensional 
30 matrix of points is obtained that defines the values that each point assumes 
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the volumetric function. 

It is possible to apply a low-pass filter to this matrix so as to obtain a 
smoother reconstructed surface removing possible errors of sampling 
without acting on the viscosity parameters (that could prejudice the correct 
5 determination of concave zones). The implementation of this filter is the 

same as that to carry out a mono and two-dimensional filtering: the entire 
volumetric function is conveyed with a three-dimensional matrix containing 
values of the filter. For example it can be made in low-pass filter using a 
matrix 3x3x3 in which the central element has the maximum value and the 

10 other elements possess a value given by the value of a Gaussian function 

depending on their distance from the centre of the matrix. 

With the present method, in which case it is possible and convenient to 
decompose a surface of large dimensions in distinct parts, it is possible to 
calculate the evolution singularly on each part and recompose then the 

1 5 various volumetric ftmctions. hi this manner it is also possible to use 

different resolutions for the different parts, then recomposing them on a 
single common matrix. It is preferable, in order to improve the link between 
the different parts, between the different defined zones, that there is 
overlapping, that is that part of the end points belonging to a zone are also 

20 incorporated in the adjacent zone. Once the matrices of the various zones are 

obtained, they are arranged on a common matrix with resolution equal to the 
maximum resolution used. In the case the contents of a cell of the final grid 
are contemporarily described by two volumetric functions (because they 
overlap) the minimum between all the volumetric functions that describe it 

25 is chosen as value for this cell. 

The method for determining the three-dimensional surface of an object, 
herein described, can be transcribed in a program code, in a manner well 
known to a technician of the sector, that can be memorized on any type of 
memory or support (floppy, CD) and/or can be carried out by a computer. 

30 



